<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Atomically resolve a metabolic reconstruction</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S4 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S11 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S12 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement,.eoOutputWrapper .matrixElement {    min-height: 18px;    box-sizing: border-box;}
.embeddedOutputsMatrixElement .matrixElement,.eoOutputWrapper  .matrixElement,.rtcDataTipElement .matrixElement {    position: relative;}
.matrixElement .variableValue,.rtcDataTipElement .matrixElement .variableValue {    white-space: normal;    display: inline-block;    vertical-align: top;    overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {    display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {    padding-top: 0 !important;    max-height: 100px;}
.inlineElement .matrixElement {    max-height: 300px;}
.embeddedOutputsMatrixElement.rightPaneElement {}
.rightPaneElement .matrixElement,.rtcDataTipElement .matrixElement {    overflow: hidden;    padding-left: 9px;}
.rightPaneElement .matrixElement {    margin-bottom: -1px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer,.eoOutputWrapper .matrixElement .valueContainer,.rtcDataTipElement .matrixElement .valueContainer {    white-space: nowrap;    margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis.hide,.eoOutputWrapper .matrixElement .verticalEllipsis.hide,.rtcDataTipElement .matrixElement .valueContainer .horizontalEllipsis.hide,.rtcDataTipElement .matrixElement .verticalEllipsis.hide {    display: none;}
.embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .verticalEllipsis, .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .horizontalEllipsis {    display:none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis {    margin-bottom: -3px;}
.eoOutputWrapper .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer {    cursor: default !important;}
/* * Ellipses as base64 for HTML export. */.matrixElement .horizontalEllipsis,.rtcDataTipElement .matrixElement .horizontalEllipsis {    display: inline-block;    margin-top: 3px;    /* base64 encoded version of images-liveeditor/HEllipsis.png */    width: 30px;    height: 12px;    background-repeat: no-repeat;    background-image: url("");}
.matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis {    margin-left: 35px;    /* base64 encoded version of images-liveeditor/VEllipsis.png */    width: 12px;    height: 30px;    background-repeat: no-repeat;    background-image: url("");}
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S16 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S17 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Atomically resolve a metabolic reconstruction</span></h1><h2  class = 'S1'><span>Author(s): </span><span style=' font-weight: bold;'>Hulda S. Haraldsdóttir and German A. Preciat Gonzalez, </span><span>Systems Biochemistry Group, University of Luxembourg.</span></h2><h2  class = 'S1'><span>Reviewer(s): Catherine Clancy, Molecular Systems Physiology Group, University of Luxembourg.</span></h2><h2  class = 'S1'><span>Francisco J. Planes, Department of Biomedical Engineering and Sciences, Tecnun, University of Navarra.</span></h2><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>Genome-scale metabolic network reconstructions have become a relevant tool in modern biology to study the metabolic pathways of biological systems </span><span style=' font-style: italic;'>in silico</span><span>. However, a more detailed representation at the underlying level of atom mappings opens the possibility for a broader range of biological, biomedical and biotechnological applications than with stoichiometry alone.</span><span> </span></div><div  class = 'S2'><span>A set of atom mappings represents the mechanism of each chemical reaction in a metabolic network, each of which relates an atom in a substrate metabolite to an atom of the same element in a product metabolite (Figure 1).</span><span> To atom map reactions in a metabolic network reconstruction, one requires chemical structures in a data file format (SMILES, MDL MOL, InChIs), reaction stoichiometries, and an atom mapping algorithm.</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "585" height = "178" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 1. Set of atom mappings for reaction L-Cysteine L-Homocysteine-Lyase (VMH ID: r0193).</span></div><div  class = 'S2'><span>Metabolites chemical structures can be obtained by different approaches such as draw them based on the literature using chemoinformatic software, or obtain them from metabolic databases either manually or using a computational software as suggested in</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;1&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. Here we recommend downloading the metabolites structures in MDL MOL format for the latest human metabolic reconstruction Recon 3</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span> via the Virtual Metabolic Human database (VMH, </span><a href = "http://vmh.life"><span>http://vmh.life</span></a><span>). Chemical structures and reaction stoichiometries from COBRA models are used to generate an MDL RXN file, which contains the information of a chemical reaction. Atom mapped reactions from Recon 3 can also be found in the VMH database in MDL RXN format. However, here we will atom map the chemical reactions using the Reaction Decoder Tool (RDT) algorithm</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;3&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>, which was selected after comparing the performance of recently published algorithms </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;4&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. However, despite its good performance (accuracy and availability) RDT algorithm does not atom map hydrogen atoms.</span></div><div  class = 'S2'><span>In this tutorial, we will identify the conserved moieties using atom mapping data for the dopamine synthesis network (DAS) extracted from Recon 3 </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span> (Figure 2). Section 1 of the tutorial will cover obtaining and visualising an atom map of metabolic reactions, and section 2 of the tutorial covers the identification of conserved metabolic moieties.</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "651" height = "658" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 2: DAS: a small metabolic network consisting of reactions in the human dopamine synthesis pathway</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. Atoms belonging to the same conserved moiety have identically coloured backgrounds.</span></div><h2  class = 'S1'><span>MATERIALS</span></h2><div  class = 'S2'><span>To atom map reactions it is required to have Java version 8 and Linux. The atom mapping does not run on Windows at present. </span></div><div  class = 'S2'><span>On </span><span style=' font-style: italic;'>macOS</span><span>, please make sure that you run the following commands in the Terminal before continuing with this tutorial:</span></div><div  class = 'S2'><span style=' font-family: monospace;'>$ /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"</span></div><div  class = 'S2'><span style=' font-family: monospace;'>$ brew install coreutils</span></div><div  class = 'S2'><span style=' font-family: monospace;'>On </span><span style=' font-style: italic; font-family: monospace;'>Linux</span><span style=' font-family: monospace;'>,</span><span>please make sure that Java and ChemAxon directories are included. To do this, run the following commands:</span></div><div  class = 'S2'><span style=' font-family: monospace;'>$ export PATH=$PATH:/opt/opt/chemaxon/jchemsuite/bin/ </span><span>(default location of JChem)</span></div><div  class = 'S2'><span style=' font-family: monospace;'>$ export PATH=$PATH:/usr/java/jre1.8.0_131/bin/ </span><span>(default installation of Java)</span></div><div  class = 'S2'><span>Also, in order to standardise the chemical reaction format it is required to have JChem downloaded from ChemAxon with its respective license.</span></div><h2  class = 'S4'><span>SECTION 1 Atom mapping of reactions</span></h2><div  class = 'S2'><span>Atom mappings for the internal reactions of a metabolic network reconstruct</span><span>ion are performed by the function </span><span style=' font-family: monospace;'>obtainAtomMappingsRDT</span><span>. The main inputs are a COBRA model structure and a directory containing the molecular structures in MDL MOL format. For this tutorial, using the RDT algorithm, the atom mappings are generated based on the molecular structures contained in cobratoolbox/tutorials/atomicallyResolveReconstruction/data/molFiles (</span><span style=' font-family: monospace;'>molFileDir</span><span>) and the reconstructed DAS network without hydrogen atoms (</span><span style=' font-family: monospace;'>model</span><span>).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">global </span><span >CBTDIR</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >tutorialdir = fileparts(which(</span><span style="color: rgb(170, 4, 249);">'tutorial_atomicallyResolveReconstruction.mlx'</span><span >));</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >model = readCbModel([tutorialdir filesep </span><span style="color: rgb(170, 4, 249);">'data' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'subDas.mat'</span><span >]); </span><span style="color: rgb(2, 128, 9);">% The subnetwork of the dopamine synthesis network</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="D61A69CC" data-testid="output_0" data-width="428" data-height="230" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">model = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">              S: [10×10 double]
           mets: {10×1 cell}
              b: [10×1 double]
         csense: [10×1 char]
           rxns: {10×1 cell}
             lb: [10×1 double]
             ub: [10×1 double]
              c: [10×1 double]
         osense: -1
          genes: {0×1 cell}
          rules: {10×1 cell}
     metCharges: [10×1 double]
    metFormulas: {10×1 cell}
    description: 'subDas.mat'
        modelID: 'model'
</div></div></div></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >molFileDir = [tutorialdir filesep </span><span style="color: rgb(170, 4, 249);">'data' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'molFiles'</span><span >]; </span><span style="color: rgb(2, 128, 9);">% The chemical structures of metabolites</span></span></div></div></div><div  class = 'S10'><span>The function </span><span style=' font-family: monospace;'>obtainAtomMappingsRDT </span><span>generates 4 different directories containing: </span></div><ul  class = 'S11'><li  class = 'S12'><span>the atom mapped reactions in MDL RXN format (directory </span><span style=' font-style: italic;'>atomMapped</span><span>), </span></li><li  class = 'S12'><span>the images of the atom mapped reactions (directory </span><span style=' font-style: italic;'>images</span><span>), </span></li><li  class = 'S12'><span>additional data for the atom mapped reactions (SMILES,  and product and reactant indexes) (directory </span><span style=' font-style: italic;'>txtData</span><span>), and </span></li><li  class = 'S12'><span>the unmapped MDL RXN files (directory </span><span style=' font-style: italic;'>rxnFiles</span><span>).</span><span style=' font-family: monospace;'> </span></li></ul><div  class = 'S2'><span>The input variable </span><span style=' font-family: monospace;'>outputDir </span><span>indicates the directory where the folders will be generated (by default the function assigns the current directory).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: normal"><span >outputDir = [pwd filesep </span><span style="color: rgb(170, 4, 249);">'output'</span><span >];</span></span></div></div></div><div  class = 'S10'><span>For some reactions, the RDT algorithm cannot compute the atom mappings (for a large reaction is generated an MDL RXN v3000 which is not compatible with the RDT algorithm). Therefore, it is necessary to assign a maximum time of processing </span><span style=' font-family: monospace;'>maxTime</span><span> (by default the function assign 30 minutes as a maximum time for computing an atom mapping for a reaction).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: normal"><span >maxTime = 1800; </span><span style="color: rgb(2, 128, 9);">% seconds</span></span></div></div></div><div  class = 'S10'><span>The function </span><span style=' font-family: monospace;'>obtainAtomMappingsRDT</span><span> generates atom mapped reactions in a standard canonical format but it is </span><span style=' font-weight: bold;'>REQUIRED</span><span> to have a Chemaxon license installed. However, the reactions can be atom mapped without being standardised. The variable </span><span style=' font-family: monospace;'>isChemaxonInstalled</span><span> contains a logical value defined by the user if the license is installed or not.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S13'><span style="white-space: normal"><span >isChemaxonInstalled = false; </span><span style="color: rgb(2, 128, 9);">% Change varibale to "true" if ChemAxon is installed</span></span></div></div></div><div  class = 'S10'><span>Now, let's obtain the atom map using </span><span style=' font-family: monospace;'>obtainAtomMappingsRDT</span><span>: </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >ispc</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    error(</span><span style="color: rgb(170, 4, 249);">'Error: atom mapping function should be run on Linux or MAC.'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    tic</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">try</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >        standardisedRxns = obtainAtomMappingsRDT(model, molFileDir, outputDir, maxTime, isChemaxonInstalled);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    toc</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="ED335C1E" data-testid="output_1" data-width="428" data-height="171" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Generating RXN files.
Computing atom mappings for 4 reactions.


4 reactions were atom mapped
4 reactions are not standardised
0 reactions were not mapped


RDT algorithm was developed by:
SA Rahman et al.: Reaction Decoder Tool (RDT): Extracting Features from Chemical
Reactions, Bioinformatics (2016), doi: 10.1093/bioinformatics/btw096</div></div><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="79BA2080" data-testid="output_2" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Elapsed time is 13.696865 seconds.</div></div></div></div></div><div  class = 'S10'><span>The output, </span><span style=' font-family: monospace;'>standardisedRxns,</span><span> is a list of atom mapped mass balanced reactions.</span></div><div  class = 'S2'><span style=' font-weight: bold;'>TIMING</span></div><div  class = 'S2'><span>The time to compute atom mappings for metabolic reactions depends on the size of the genome-sca</span><span>le model and the size of the molecules in the reactions. The above example may take ~1 min or less if </span><span style=' font-family: monospace;'>isChemaxonInstalled = false.</span></div><h2  class = 'S4'><span>Visualising results</span></h2><div  class = 'S2'><span>The </span><span style=' font-style: italic;'>images</span><span> directory contains a graphical representation of the atom mapped reactions. They show the bijection between atoms and each of the metabolite pools are coloured for an easy visualisation. Figure 3 shows the atom mapped reaction to produce dopamine and </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span> from L-DOPA.</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "210" height = "448" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 3: Reaction 3-hydroxy-L-tyrosine carboxy-lyase atom mapped (VMH ID: 3HLYCL) here represented as R3. Images ge</span><span>nerated by RDT algorithm, also shows where a reaction centre occurs.</span></div><div  class = 'S2'><span>The </span><span style=' font-style: italic;'>rxnFiles </span><span>directory contains for all atom mapped reactions a corresponding MDL RXN file (Figure 4). Contained within these files are information of the chemical reaction, such as: </span></div><ul  class = 'S11'><li  class = 'S12'><span>the name of the reaction (on line 2 of the file), </span></li><li  class = 'S12'><span>the chemical formula (on line 4 of the file), </span></li><li  class = 'S12'><span>the number of substrates and products (on line 5 of the file), and </span></li><li  class = 'S12'><span>specific information for each of the molecules (from line 6 onwards, after the identifier $MOL). </span></li></ul><div  class = 'S3'><img class = "imageNode" src = "" width = "505" height = "545" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 4: A MDL RXN file stored in the </span><span style=' font-style: italic;'>rxnFiles </span><span>directory.  </span></div><div  class = 'S2'><span>Specific information for each of the molecules includes the name of the metabolite, it's InChI key (if the metabolite does not contain an R group) and the number of atoms and bonds. Following this is the atom block, which contains detailed information on the coordinates, element, charge and atom mapping number for each of the atoms, and then finally, the bond block connects all the atoms in the metabolite.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >regexp(fileread([outputDir filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'R3.rxn'</span><span >]), </span><span style="color: rgb(170, 4, 249);">'\n'</span><span >, </span><span style="color: rgb(170, 4, 249);">'split'</span><span >)'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="6EA920F4" data-testid="output_3" data-width="428" data-height="1168" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    '$RXN'
    ''
    '  EC-BLAST     R3'
    ''
    '  1  2'
    '$MOL'
    '34dhphe[c]'
    '  EC-BLAST  1113171051'
    'InChIKey=WTDRDQBEARUVNC-LURJTMIESA-N'
    ' 14 14  0  0  0  0  0  0  0  0999 V2000'
    '    0.7145    3.7125    0.0000 O   0  0  0  0  0  0  0  0  0  1  0  0'
    '    0.0000    3.3000    0.0000 C   0  0  0  0  0  0  0  0  0  2  0  0'
    '   -0.7145    3.7125    0.0000 O   0  0  0  0  0  0  0  0  0  3  0  0'
    '    0.0000    2.4750    0.0000 C   0  0  1  0  0  0  0  0  0  4  0  0'
    '   -0.7145    2.0625    0.0000 N   0  0  0  0  0  0  0  0  0  5  0  0'
    '    0.7145    2.0625    0.0000 C   0  0  0  0  0  0  0  0  0  6  0  0'
    '    0.7145    1.2375    0.0000 C   0  0  0  0  0  0  0  0  0  7  0  0'
    '    1.4289    0.8250    0.0000 C   0  0  0  0  0  0  0  0  0  8  0  0'
    '    1.4289    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  9  0  0'
    '    0.7145   -0.4125    0.0000 C   0  0  0  0  0  0  0  0  0 10  0  0'
    '    0.7145   -1.2375    0.0000 O   0  0  0  0  0  0  0  0  0 11  0  0'
    '    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0 12  0  0'
    '   -0.7145   -0.4125    0.0000 O   0  0  0  0  0  0  0  0  0 13  0  0'
    '    0.0000    0.8250    0.0000 C   0  0  0  0  0  0  0  0  0 14  0  0'
    '  2  1  2  0  0  0  0 '
    '  2  3  1  0  0  0  0 '
    '  4  2  1  0  0  0  0 '
    '  4  5  1  6  0  0  0 '
    '  4  6  1  0  0  0  0 '
    '  6  7  1  0  0  0  0 '
    '  7  8  2  0  0  0  0 '
    '  7 14  1  0  0  0  0 '
    '  8  9  1  0  0  0  0 '
    '  9 10  2  0  0  0  0 '
    ' 10 11  1  0  0  0  0 '
    ' 10 12  1  0  0  0  0 '
    ' 12 13  1  0  0  0  0 '
    ' 12 14  2  0  0  0  0 '
    'M  CHG  1   3  -1'
    'M  CHG  1   5   1'
    'M  END'
    '$MOL'
    'dopa[c]'
    '  EC-BLAST  1113171051'
    'InChIKey=VYFYYTLLBUKUHU-UHFFFAOYSA-O'
    ' 11 11  0  0  0  0  0  0  0  0999 V2000'
    '   -0.7145   -0.4125    0.0000 O   0  0  0  0  0  0  0  0  0 13  0  0'
    '    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0 12  0  0'
    '    0.0000    0.8250    0.0000 C   0  0  0  0  0  0  0  0  0 14  0  0'
    '    0.7145    1.2375    0.0000 C   0  0  0  0  0  0  0  0  0  7  0  0'
    '    1.4289    0.8250    0.0000 C   0  0  0  0  0  0  0  0  0  8  0  0'
    '    0.0000    3.3000    0.0000 N   0  0  0  0  0  0  0  0  0  5  0  0'
    '    0.7145    2.0625    0.0000 C   0  0  0  0  0  0  0  0  0  6  0  0'
    '    0.0000    2.4750    0.0000 C   0  0  0  0  0  0  0  0  0  4  0  0'
    '    1.4289    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  9  0  0'
    '    0.7145   -1.2375    0.0000 O   0  0  0  0  0  0  0  0  0 11  0  0'
    '    0.7145   -0.4125    0.0000 C   0  0  0  0  0  0  0  0  0 10  0  0'
    '  2  1  1  0  0  0  0 '
    '  2  3  2  0  0  0  0 '
    ' 11  2  1  0  0  0  0 '
    '  4  3  1  0  0  0  0 '
    '  4  5  2  0  0  0  0 '
    '  7  4  1  0  0  0  0 '
    '  5  9  1  0  0  0  0 '
    '  6  8  1  0  0  0  0 '
    '  8  7  1  0  0  0  0 '
    '  9 11  2  0  0  0  0 '
    ' 11 10  1  0  0  0  0 '
    'M  CHG  1   6   1'
    'M  END'
    '$MOL'
    'co2[c]'
    '  EC-BLAST  1113171051'
    'InChIKey=CURLTUGMZLYLDI-UHFFFAOYSA-N'
    '  3  2  0  0  0  0  0  0  0  0999 V2000'
    '    1.6500    0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  1  0  0'
    '    2.4750    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  2  0  0'
    '    3.3000    0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  3  0  0'
    '  1  2  2  0  0  0  0 '
    '  2  3  2  0  0  0  0 '
    'M  END'
    ''
</div></div></div></div></div></div><div  class = 'S10'><span>The </span><span style=' font-style: italic;'>txtData </span><span>directory contains the TXT information of the reaction including the SMILES format, which holds the standard canonical format of the reaction, the reactant input atom index and the product input atom index.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >regexp(fileread([outputDir filesep </span><span style="color: rgb(170, 4, 249);">'txtData' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'R3.txt'</span><span >]), </span><span style="color: rgb(170, 4, 249);">'\n'</span><span >, </span><span style="color: rgb(170, 4, 249);">'split'</span><span >)'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="9C6EF6EB" data-testid="output_4" data-width="428" data-height="202" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    ''
    '//'
    'SELECTED AAM MAPPING'
    '[O:1]=[C:2]([O-:3])[CH:4]([NH3+:5])[CH2:6][C:7]:1:[CH:8]:[CH:9]:[C:10]([OH:11]):[C:12]([OH:13]):[CH:14]1&gt;&gt;[O:1]=[C:2]=[O:3].[OH:11][C:10]:1:[CH:9]:[CH:8]:[C:7](:[CH:14]:[C:12]1[OH:13])[CH2:6][CH2:4][NH3+:5]'
    ''
    ''
    '//'
    'REACTANT INPUT ATOM INDEX&lt;-&gt;AAM ID'
    '{1=5, 2=4, 3=6, 4=7, 5=8, 6=9, 7=10, 8=11, 9=12, 10=13, 11=14, 12=2, 13=3, 14=1}'
    'PRODUCT INPUT ATOM INDEX&lt;-&gt;AAM ID'
    '{1=11, 2=10, 3=9, 4=5, 5=4, 6=3, 7=2, 8=1, 9=7, 10=8, 11=6, 12=12, 13=13, 14=14}'
    ''
    ''
</div></div></div></div></div></div><h2  class = 'S4'><span>SECTION 2 Identifying conserved metabolic moieties</span></h2><div  class = 'S2'><span>A conserved moiety is a group of atoms within molecules connected by reactions, that follow identical paths through a metabolic network and therefore, its amount remains constant (Figure 5). Representative examples from energy metabolism include the AMP and NAD moieties. With the set of atom mappings for a metabolic network the set of linearly independent conserved moieties for the metabolic network can be identified, each of which corresponds to a particular identifiable molecular substructure</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;5&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>.</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "425" height = "194" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 5: A graphical representation of a conserved moiety</span></div><div  class = 'S2'><span>In this section, we will identify conserved moieties in a subnetwork of the DAS network (Figure 2) by graph theoretical analysis of its atom transition network. The method is described in</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;5&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>. This section consists of two parts: </span></div><div  class = 'S2'><span>Part 1 covers basic usage of the code. </span></div><div  class = 'S2'><span>Part 2 covers decomposition of a composite moiety resulting from variable atom mappings between the recurring metabolite pair </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;O&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="19" height="20" /></span><span> and </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;H&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;O&lt;/mi&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="29.5" height="20" /></span><span>.</span></div><h2  class = 'S4'><span>Part 1: Identify conserved moieties in DAS</span></h2><div  class = 'S2'><span style=' font-weight: bold;'>Step 1: Generate an atom transition network for DAS based on atom mappings for internal (mass and charge balanced) reactions.</span></div><div  class = 'S2'><span>The atom transition network is generated based on the reconstructed DAS network (</span><span style=' font-family: monospace;'>model</span><span>) and atom mappings for internal reactions, obtained in the previous section and predicted with the RDT algorithm</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;3&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~isChemaxonInstalled</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    copyfile([tutorialdir filesep </span><span style="color: rgb(170, 4, 249);">'data' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped'</span><span >],[outputDir filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped'</span><span >])</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span><span >    </span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >atomMappedDir = [outputDir filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped'</span><span >];</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >ATN = buildAtomTransitionNetwork(model, atomMappedDir);</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="60F0F92A" data-testid="output_5" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Atom mappings found for 4 model reactions.
Generating atom transition network for reactions with atom mappings.</div></div></div></div></div><div  class = 'S10'><span>The output variable (</span><span style=' font-family: monospace;'>ATN</span><span>) is a Matlab structure with several fields. </span><span style=' font-family: monospace;'>ATN.A</span><span> is </span><span>the incidence matrix of the directed graph representing the atom transition network. Each row represents a particular atom in one of the 11 DAS metabolites. </span><span style=' font-family: monospace;'>ATN.mets</span><span> indicates which metabolite in DAS each atom belongs to. To find rows of </span><span style=' font-family: monospace;'>ATN.A </span><span>corresponding to atoms in </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span>, run:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >ico2 = find(ismember(ATN.mets, </span><span style="color: rgb(170, 4, 249);">'co2[c]'</span><span >))'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="8DD3F3D2" data-testid="output_6" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>ico2 = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:43.2,&quot;totalColumns&quot;:&quot;3&quot;,&quot;totalRows&quot;:&quot;1&quot;,&quot;charsPerColumn&quot;:6}"><div class="variableValue" style="width: 131.6px;">    88    89    90
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><div  class = 'S10'><span>The order of atoms in </span><span style=' font-family: monospace;'>ATN.A</span><span> matches their order in MDL MOL files encoding metabolite structures (Figure 7), e.g., </span><span style=' font-family: monospace;'>ATN.A</span><span>(90,:) is the row corresponding to the second oxygen atom (number 3 in Figure 6).</span></div><div  class = 'S2'><span></span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "114" height = "33" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 6: Rows for </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span> atoms in </span><span style=' font-family: monospace;'>ATN.A</span><span> are ordered as shown.</span></div><div  class = 'S2'><span style=' font-family: monospace;'></span></div><div  class = 'S2'><span style=' font-family: monospace;'>ATN.elements</span><span> contains the element symbols of atoms, e.g.,</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >ATN.elements{90}</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="1E58D9A7" data-testid="output_7" data-width="428" data-height="20" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span>O</div></div></div></div></div></div><div  class = 'S10'><span>Each column of </span><span style=' font-family: monospace;'>ATN.A</span><span> represents a particular atom transition in one of the four internal reactions in DAS. Reaction identifiers of atom transitions are given in </span><span style=' font-family: monospace;'>ATN.rxns</span><span>. To find all atom transitions that involve </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span> atoms, run:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >tco2 = find(any(ATN.A(ico2,:), 1))</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="82ADB73D" data-testid="output_8" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>tco2 = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:43.2,&quot;totalColumns&quot;:&quot;6&quot;,&quot;totalRows&quot;:&quot;1&quot;,&quot;charsPerColumn&quot;:6}"><div class="variableValue" style="width: 261.2px;">    75    76    77    95    96    97
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >ATN.rxns(tco2)'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="CCB76B9D" data-testid="output_9" data-width="428" data-height="34" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'R3'    'R3'    'R3'    'R4'    'R4'    'R4'
</div></div></div></div></div></div><div  class = 'S10'><span>i.e., three atom transitions in each of the reactions R3 and R4 involve atoms in </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span>. To find atoms connected to </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span> atoms via these atom transitions, run:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >cco2 = find(any(ATN.A(:, tco2) &lt; 0,2));</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >ATN.mets(cco2)'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="DD2736C6" data-testid="output_10" data-width="428" data-height="34" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    '34dhphe[c]'    '34dhphe[c]'    '34dhphe[c]'    'for[c]'    'for[c]'    'for[c]'
</div></div></div></div></div></div><div  class = 'S10'><span>i.e., </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot;&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;CO&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;2&lt;/mn&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="28.5" height="20" /></span><span> atoms are connected to atoms in the metabolites L-DOPA (VMH ID: 34dhphe) and formate (VMH ID: for).</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Step 2: Identify conserved moieties in DAS by graph theoretical analysis of the atom transition network generated in Step 1.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >tic</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >[L,Lambda,moietyFormulas,moieties2mets,moieties2vectors,atoms2moieties] = </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    identifyConservedMoieties(model, ATN);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >t = toc;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >fprintf(</span><span style="color: rgb(170, 4, 249);">'Computation time: %.1e s\n\n'</span><span >, t); </span><span style="color: rgb(2, 128, 9);">% Print computation time</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="6E7A2230" data-testid="output_11" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Computation time: 3.2e-01 s</div></div></div></div></div><div  class = 'S10'><span>This function outputs the moiety matrix (</span><span style=' font-family: monospace;'>L</span><span>), the moiety supergraph (</span><span style=' font-family: monospace;'>Lambda</span><span>), the chemical formulas of moieties (</span><span style=' font-family: monospace;'>moietyFormulas</span><span>), and three vectors that map between the various inputs and outputs. The 10×5 moiety matrix </span><span style=' font-family: monospace;'>L</span><span> has a row for each metabolite and a column for each conserved moiety in DAS. Each column is a moiety vector, with elements corresponding to the number of instances of a conserved moiety in each metabolite. To find the number of instances of moiety 2 in L-DOPA, run</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >iLDOPA = find(ismember(model.mets, </span><span style="color: rgb(170, 4, 249);">'34dhphe[c]'</span><span >))</span></span></div><div  class = 'S8'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>iLDOPA = 7</div></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >full(L(iLDOPA, 2))</span></span></div><div  class = 'S8'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 1</div></div></div></div><div  class = 'S10'><span>i.e., L-DOPA contains one instance of moiety 2.</span></div><div  class = 'S2'><span>The 19×17 moiety supergraph (</span><span style=' font-family: monospace;'>Lambda</span><span>) contains the graphs of all seven conserved moieties in DAS (Figure 7).</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "517" height = "285" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 7: Graphs of the five conserved moieties in DAS. Each node represents an instance of a conserved moiety in a particular metabolite. Each directed edge represents conservation of a moiety between two metabolites. The chemical formula of each moiety is given below its graph.</span></div><div  class = 'S2'><span>Each row of Lambda represents a single instance of a conserved moiety in a particular metabolite. The vector moieties2vectors maps between the rows of Lambda and the columns of L. To obtain the incidence matrix of a particular moiety graph, e.g., λ2 in Figure 7, run</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >i2 = find(moieties2vectors == 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >c2 = find(any(Lambda(i2, :)));</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >lambda2 = full(Lambda(i2, c2))</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="EDDB1083" data-testid="output_14" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>lambda2 = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:43.2,&quot;totalColumns&quot;:&quot;4&quot;,&quot;totalRows&quot;:&quot;5&quot;,&quot;charsPerColumn&quot;:6}"><div class="variableValue" style="width: 174.8px;">    -1     0     0     0
     1    -1     0     0
     0     1    -1     0
     0     0     1     1
     0     0     0    -1
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><div  class = 'S10'><span>The vector </span><span style=' font-family: monospace;'>moieties2mets</span><span> maps the rows of Lambda to metabolite indices in the DAS reconstruction (</span><span style=' font-family: monospace;'>model</span><span>). To find metabolites containing instances of moiety 2, run</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >m2 = moieties2mets(i2);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >mets2 = model.mets(m2)'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="81E19FB8" data-testid="output_15" data-width="428" data-height="34" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">mets2 = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'phe_L[c]'    'tyr_L[c]'    '34dhphe[c]'    'co2[c]'    'for[c]'
</div></div></div></div></div></div><div  class = 'S10'><span>The chemical formula of moiety 2 is given by,</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >moietyFormulas{2}</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="6EDD72C1" data-testid="output_16" data-width="428" data-height="20" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span>CO2</div></div></div></div></div></div><div  class = 'S10'><span>Finally, the vector atoms2moieties maps each atom in the atom transition network for DAS to a particular instance of a conserved moiety. To find atoms in L-DOPA that belong to moiety 2, run</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >find(ismember(atoms2moieties, i2) &amp; ismember(ATN.mets, </span><span style="color: rgb(170, 4, 249);">'34dhphe[c]'</span><span >))'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="DABB63DC" data-testid="output_17" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>ans = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:43.2,&quot;totalColumns&quot;:&quot;3&quot;,&quot;totalRows&quot;:&quot;1&quot;,&quot;charsPerColumn&quot;:6}"><div class="variableValue" style="width: 131.6px;">    74    75    76
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><div  class = 'S10'><span style=' font-weight: bold;'>Step 3: Classify moieties</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >types = classifyMoieties(L, model.S)</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="C1DF56D0" data-testid="output_18" data-width="428" data-height="90" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">types = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'Transitive'
    'Transitive'
    'Internal'
    'Transitive'
    'Transitive'
</div></div></div></div></div></div><div  class = 'S10'><span>The internal moiety (λ3 in Figure 3) is conserved in both the open and closed DAS network, whereas the transitive and integrative moieties are only conserved in the closed network</span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mn&gt;6&lt;/mn&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="11" height="19" /></span><span>.</span></div><h2  class = 'S4'><span>Part 2: Effects of variable atom mappings between recurring metabolite pairs</span></h2><div  class = 'S2'><span>Here, we will again identify conserved moieties in DAS but with a slightly different set of atom mappings (Figure 8). The different atom mappings gives rise to a different atom transition network with a different set of conserved moieties. In particular, it contains a single composite moiety, λ8 in Figure 5, in place of the two moieties λ4 and λ5 in Figure 3. The composite moiety is the result of variable atom mappings between the recurring metabolite pair O2 and H2O in reactions R1 and R2.</span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "549" height = "260" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>Figure 8: (a) Oxygen atom transitions used in Part 1. Oxygen atom 1 in O2 maps to the oxygen atom in H2O in both R1 and R2. These atom transitions contain two separate moieties, with two disconnected moiety graphs (λ4 and λ5 in Figure 7), and two linearly independent moiety vectors (L(:,4) and L(:,5)). (b) Oxygen atom transitions used in Part 2. A different oxygen atom maps from O2 to H2O in R1 than in R2. These atom transitions contain only one composite moiety. (c) The composite moiety graph arising from the oxygen atom transitions in (b).</span></div><div  class = 'S2'><span style=' font-weight: bold;'>Step 1: Identify conserved moieties with the alternative set of atom mappings.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% Create an alternative MDL RXN file</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >R2rxn = regexp(fileread([outputDir filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'R2.rxn'</span><span >]), </span><span style="color: rgb(170, 4, 249);">'\n'</span><span >, </span><span style="color: rgb(170, 4, 249);">'split'</span><span >)';</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >R2rxn{2} = </span><span style="color: rgb(170, 4, 249);">'alternativeR2'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >R2rxn{135}(62:63) = </span><span style="color: rgb(170, 4, 249);">'18'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >R2rxn{151}(62:63) = </span><span style="color: rgb(170, 4, 249);">'19'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >fid2 = fopen([outputDir filesep </span><span style="color: rgb(170, 4, 249);">'atomMapped' </span><span >filesep </span><span style="color: rgb(170, 4, 249);">'alternativeR2.rxn'</span><span >], </span><span style="color: rgb(170, 4, 249);">'w'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >fprintf(fid2, </span><span style="color: rgb(170, 4, 249);">'%s\n'</span><span >, R2rxn{:});</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >fclose(fid2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% Create an alternative DAS model</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >alternativeModel = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >alternativeModel.rxns{2} = </span><span style="color: rgb(170, 4, 249);">'alternativeR2'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% Identify conserved moieties</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >ATN = buildAtomTransitionNetwork(alternativeModel, atomMappedDir);</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="DA1A124A" data-testid="output_19" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">Atom mappings found for 4 model reactions.
Generating atom transition network for reactions with atom mappings.</div></div></div></div><div class="inlineWrapper"><div  class = 'S16'><span style="white-space: normal"><span >[L,Lambda,moietyFormulas,moieties2mets,moieties2vectors,atoms2moieties] = </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S17'><span style="white-space: normal"><span >    identifyConservedMoieties(alternativeModel, ATN);</span></span></div></div></div><div  class = 'S2'><span style=' font-weight: bold;'>Step 2: Decompose the composite moiety vector</span></div><div  class = 'S2'><span>First, extract the internal stoichiometric matrix for DAS, by running:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >rbool = ismember(alternativeModel.rxns, ATN.rxns);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >mbool = any(alternativeModel.S(:,rbool), 2);</span></span></div></div><div class="inlineWrapper"><div  class = 'S17'><span style="white-space: normal"><span >N = alternativeModel.S(mbool, rbool);</span></span></div></div></div><div  class = 'S10'><span>To decompose the moiety matrix computed in Step 1, run:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">try</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'gurobi6'</span><span >, </span><span style="color: rgb(170, 4, 249);">'milp'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="5F18BAC0" data-testid="output_20" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; Gurobi interface added to MATLAB path.
 &gt; gurobi (version 751) is compatible and fully tested with MATLAB R2016b on your operating system.</div></div></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >D = decomposeMoietyVectors(L, N);</span></span></div></div></div><div  class = 'S10'><span>Note that you can use any Mixed Integer Linear Programme (MILP) solver that is supported by the COBRA toolbox. The decomposed moiety matrix D is identical to the original moiety matrix computed in Part 1. Moiety vectors D(:,4) and D(:,5) are the linearly independent components of the composite moiety vector L(:,4) above.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >full(D(:,[4 5])')</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="E943B889" data-testid="output_21" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>ans = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:43.2,&quot;totalColumns&quot;:&quot;10&quot;,&quot;totalRows&quot;:&quot;2&quot;,&quot;charsPerColumn&quot;:6}"><div class="variableValue" style="width: 390.8px;">     0     0     1     0     0     1     0     0     0     0
     0     0     1     1     0     0     2     2     0     0
</div><div class="horizontalEllipsis"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><div  class = 'S10'><span>One disadvantage of decomposing moiety vectors is that it is difficult to keep track of which atoms belong to the decomposed moieties. We can, however, estimate the chemical formulas of the decomposed moieties using the elemental matrix for DAS. The elemental matrix is a numerical representation of the chemical formulas of metabolites in DAS.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >[E,elements] = constructElementalMatrix(alternativeModel.metFormulas,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >    alternativeModel.metCharges);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: normal"><span >decomposedMoietyFormulas = estimateMoietyFormulas(D, E, elements);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S7'><span style="white-space: normal"><span >decomposedMoietyFormulas([4 5])'</span></span></div><div  class = 'S8'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="4171FFF3" data-testid="output_22" data-width="428" data-height="34" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ans = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'HO'    'O'
</div></div></div></div></div></div><div  class = 'S10'><span>i.e., each decomposed moiety contains an oxygen atom.</span></div><h2  class = 'S4'><span>References</span></h2><ol  class = 'S11'><li  class = 'S12'><span>Haraldsdóttir, H.S., Thiele, I., Fleming, R.M. Comparative evaluation of open source software for mapping between metabolite identifiers in metabolic network reconstructions: application to Recon 2. </span><span style=' font-style: italic;'>J. Cheminform</span><span> 6(1), 2 (2014).</span></li><li  class = 'S12'><span>Elizabeth Brunk, et al.</span><span> Recon 3D: A Three-Dimensional View of Human Metabolism and Disease. Submited</span></li><li  class = 'S12'><span>Rahman, S.A., et al. Reaction Decoder Tool (RDT): extracting features from chemical reactions.  </span><span style=' font-style: italic;'>Bioinformatics</span><span> 32(13), 2065–2066 (2016).</span></li><li  class = 'S12'><span>Preciat et al. Comparative evaluation of atom mapping algorithms for balanced metabolic reactions: application to Recon 3D. </span><span style=' font-style: italic;'>J Cheminform</span><span>, 9: 39 (2017).</span></li><li  class = 'S12'><span>Hulda S. Haraldsdóttir and Ronan M. T. Fleming. Identification of conserved moieties in metabolic networks by graph theoretical analysis of atom transition networks. </span><span style=' font-style: italic;'>PLOS Comput. Biol</span><span>, 12(11) (2016).</span></li><li  class = 'S12'><span>Iman Famili and B. Ø. Palsson. The convex basis of the left null space of the stoichiometric matrix leads to the definition of metabolically meaningful pools. ‎</span><span style=' font-style: italic;'>Biophys. J</span><span>, 85(1):16–26 (2003).</span></li></ol>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Atomically resolve a metabolic reconstruction
%% Author(s): *Hulda S. Haraldsdóttir and German A. Preciat Gonzalez,* Systems Biochemistry Group, University of Luxembourg.
%% Reviewer(s): Catherine Clancy, Molecular Systems Physiology Group, University of Luxembourg.
%% Francisco J. Planes, Department of Biomedical Engineering and Sciences, Tecnun, University of Navarra.
%% INTRODUCTION
% Genome-scale metabolic network reconstructions have become a relevant tool 
% in modern biology to study the metabolic pathways of biological systems _in 
% silico_. However, a more detailed representation at the underlying level of 
% atom mappings opens the possibility for a broader range of biological, biomedical 
% and biotechnological applications than with stoichiometry alone. 
% 
% A set of atom mappings represents the mechanism of each chemical reaction 
% in a metabolic network, each of which relates an atom in a substrate metabolite 
% to an atom of the same element in a product metabolite (Figure 1). To atom map 
% reactions in a metabolic network reconstruction, one requires chemical structures 
% in a data file format (SMILES, MDL MOL, InChIs), reaction stoichiometries, and 
% an atom mapping algorithm.
% 
% 
% 
% Figure 1. Set of atom mappings for reaction L-Cysteine L-Homocysteine-Lyase 
% (VMH ID: r0193).
% 
% Metabolites chemical structures can be obtained by different approaches such 
% as draw them based on the literature using chemoinformatic software, or obtain 
% them from metabolic databases either manually or using a computational software 
% as suggested in${\;}^1$. Here we recommend downloading the metabolites structures 
% in MDL MOL format for the latest human metabolic reconstruction Recon 3${\;}^2$ 
% via the Virtual Metabolic Human database (VMH, <http://vmh.life http://vmh.life>). 
% Chemical structures and reaction stoichiometries from COBRA models are used 
% to generate an MDL RXN file, which contains the information of a chemical reaction. 
% Atom mapped reactions from Recon 3 can also be found in the VMH database in 
% MDL RXN format. However, here we will atom map the chemical reactions using 
% the Reaction Decoder Tool (RDT) algorithm${\;}^3$, which was selected after 
% comparing the performance of recently published algorithms ${\;}^4$. However, 
% despite its good performance (accuracy and availability) RDT algorithm does 
% not atom map hydrogen atoms.
% 
% In this tutorial, we will identify the conserved moieties using atom mapping 
% data for the dopamine synthesis network (DAS) extracted from Recon 3 ${\;}^2$ 
% (Figure 2). Section 1 of the tutorial will cover obtaining and visualising an 
% atom map of metabolic reactions, and section 2 of the tutorial covers the identification 
% of conserved metabolic moieties.
% 
% 
% 
% Figure 2: DAS: a small metabolic network consisting of reactions in the human 
% dopamine synthesis pathway${\;}^2$. Atoms belonging to the same conserved moiety 
% have identically coloured backgrounds.
%% MATERIALS
% To atom map reactions it is required to have Java version 8 and Linux. The 
% atom mapping does not run on Windows at present. 
% 
% On _macOS_, please make sure that you run the following commands in the Terminal 
% before continuing with this tutorial:
% 
% |$ /usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"|
% 
% |$ brew install coreutils|
% 
% |On _Linux_,|please make sure that Java and ChemAxon directories are included. 
% To do this, run the following commands:
% 
% |$ export PATH=$PATH:/opt/opt/chemaxon/jchemsuite/bin/| (default location 
% of JChem)
% 
% |$ export PATH=$PATH:/usr/java/jre1.8.0_131/bin/| (default installation of 
% Java)
% 
% Also, in order to standardise the chemical reaction format it is required 
% to have JChem downloaded from ChemAxon with its respective license.
%% SECTION 1 Atom mapping of reactions
% Atom mappings for the internal reactions of a metabolic network reconstruction 
% are performed by the function |obtainAtomMappingsRDT|. The main inputs are a 
% COBRA model structure and a directory containing the molecular structures in 
% MDL MOL format. For this tutorial, using the RDT algorithm, the atom mappings 
% are generated based on the molecular structures contained in cobratoolbox/tutorials/atomicallyResolveReconstruction/data/molFiles 
% (|molFileDir|) and the reconstructed DAS network without hydrogen atoms (|model|).

global CBTDIR
tutorialdir = fileparts(which('tutorial_atomicallyResolveReconstruction.mlx'));
model = readCbModel([tutorialdir filesep 'data' filesep 'subDas.mat']); % The subnetwork of the dopamine synthesis network
molFileDir = [tutorialdir filesep 'data' filesep 'molFiles']; % The chemical structures of metabolites
%% 
% The function |obtainAtomMappingsRDT| generates 4 different directories containing: 
%% 
% * the atom mapped reactions in MDL RXN format (directory _atomMapped_), 
% * the images of the atom mapped reactions (directory _images_), 
% * additional data for the atom mapped reactions (SMILES,  and product and 
% reactant indexes) (directory _txtData_), and 
% * the unmapped MDL RXN files (directory _rxnFiles_). 
%% 
% The input variable |outputDir| indicates the directory where the folders will 
% be generated (by default the function assigns the current directory).

outputDir = [pwd filesep 'output'];
%% 
% For some reactions, the RDT algorithm cannot compute the atom mappings (for 
% a large reaction is generated an MDL RXN v3000 which is not compatible with 
% the RDT algorithm). Therefore, it is necessary to assign a maximum time of processing 
% |maxTime| (by default the function assign 30 minutes as a maximum time for computing 
% an atom mapping for a reaction).

maxTime = 1800; % seconds
%% 
% The function |obtainAtomMappingsRDT| generates atom mapped reactions in a 
% standard canonical format but it is *REQUIRED* to have a Chemaxon license installed. 
% However, the reactions can be atom mapped without being standardised. The variable 
% |isChemaxonInstalled| contains a logical value defined by the user if the license 
% is installed or not.

isChemaxonInstalled = false; % Change varibale to "true" if ChemAxon is installed
%% 
% Now, let's obtain the atom map using |obtainAtomMappingsRDT|: 

if ispc
    error('Error: atom mapping function should be run on Linux or MAC.')
else
    tic
    try
        standardisedRxns = obtainAtomMappingsRDT(model, molFileDir, outputDir, maxTime, isChemaxonInstalled);
    end
    toc
end
%% 
% The output, |standardisedRxns,| is a list of atom mapped mass balanced reactions.
% 
% *TIMING*
% 
% The time to compute atom mappings for metabolic reactions depends on the size 
% of the genome-scale model and the size of the molecules in the reactions. The 
% above example may take ~1 min or less if |isChemaxonInstalled = false.|
%% Visualising results
% The _images_ directory contains a graphical representation of the atom mapped 
% reactions. They show the bijection between atoms and each of the metabolite 
% pools are coloured for an easy visualisation. Figure 3 shows the atom mapped 
% reaction to produce dopamine and ${\textrm{CO}}_2$ from L-DOPA.
% 
% 
% 
% Figure 3: Reaction 3-hydroxy-L-tyrosine carboxy-lyase atom mapped (VMH ID: 
% 3HLYCL) here represented as R3. Images generated by RDT algorithm, also shows 
% where a reaction centre occurs.
% 
% The _rxnFiles_ directory contains for all atom mapped reactions a corresponding 
% MDL RXN file (Figure 4). Contained within these files are information of the 
% chemical reaction, such as: 
%% 
% * the name of the reaction (on line 2 of the file), 
% * the chemical formula (on line 4 of the file), 
% * the number of substrates and products (on line 5 of the file), and 
% * specific information for each of the molecules (from line 6 onwards, after 
% the identifier $MOL). 
%% 
% 
% 
% Figure 4: A MDL RXN file stored in the _rxnFiles_ directory.  
% 
% Specific information for each of the molecules includes the name of the metabolite, 
% it's InChI key (if the metabolite does not contain an R group) and the number 
% of atoms and bonds. Following this is the atom block, which contains detailed 
% information on the coordinates, element, charge and atom mapping number for 
% each of the atoms, and then finally, the bond block connects all the atoms in 
% the metabolite.

regexp(fileread([outputDir filesep 'atomMapped' filesep 'R3.rxn']), '\n', 'split')'
%% 
% The _txtData_ directory contains the TXT information of the reaction including 
% the SMILES format, which holds the standard canonical format of the reaction, 
% the reactant input atom index and the product input atom index.

regexp(fileread([outputDir filesep 'txtData' filesep 'R3.txt']), '\n', 'split')'
%% SECTION 2 Identifying conserved metabolic moieties
% A conserved moiety is a group of atoms within molecules connected by reactions, 
% that follow identical paths through a metabolic network and therefore, its amount 
% remains constant (Figure 5). Representative examples from energy metabolism 
% include the AMP and NAD moieties. With the set of atom mappings for a metabolic 
% network the set of linearly independent conserved moieties for the metabolic 
% network can be identified, each of which corresponds to a particular identifiable 
% molecular substructure${\;}^5$.
% 
% 
% 
% Figure 5: A graphical representation of a conserved moiety
% 
% In this section, we will identify conserved moieties in a subnetwork of the 
% DAS network (Figure 2) by graph theoretical analysis of its atom transition 
% network. The method is described in${\;}^5$. This section consists of two parts: 
% 
% Part 1 covers basic usage of the code. 
% 
% Part 2 covers decomposition of a composite moiety resulting from variable 
% atom mappings between the recurring metabolite pair $O_2$ and $H_2 O$.
%% Part 1: Identify conserved moieties in DAS
% *Step 1: Generate an atom transition network for DAS based on atom mappings 
% for internal (mass and charge balanced) reactions.*
% 
% The atom transition network is generated based on the reconstructed DAS network 
% (|model|) and atom mappings for internal reactions, obtained in the previous 
% section and predicted with the RDT algorithm${\;}^3$.

if ~isChemaxonInstalled
    copyfile([tutorialdir filesep 'data' filesep 'atomMapped'],[outputDir filesep 'atomMapped'])
end    
atomMappedDir = [outputDir filesep 'atomMapped'];
ATN = buildAtomTransitionNetwork(model, atomMappedDir);
%% 
% The output variable (|ATN|) is a Matlab structure with several fields. |ATN.A| 
% is the incidence matrix of the directed graph representing the atom transition 
% network. Each row represents a particular atom in one of the 11 DAS metabolites. 
% |ATN.mets| indicates which metabolite in DAS each atom belongs to. To find rows 
% of |ATN.A| corresponding to atoms in ${\textrm{CO}}_2$, run:

ico2 = find(ismember(ATN.mets, 'co2[c]'))'
%% 
% The order of atoms in |ATN.A| matches their order in MDL MOL files encoding 
% metabolite structures (Figure 7), e.g., |ATN.A|(90,:) is the row corresponding 
% to the second oxygen atom (number 3 in Figure 6).
% 
% 
% 
% 
% 
% Figure 6: Rows for ${\textrm{CO}}_2$ atoms in |ATN.A| are ordered as shown.
% 
% 
% 
% |ATN.elements| contains the element symbols of atoms, e.g.,

ATN.elements{90}
%% 
% Each column of |ATN.A| represents a particular atom transition in one of the 
% four internal reactions in DAS. Reaction identifiers of atom transitions are 
% given in |ATN.rxns|. To find all atom transitions that involve ${\textrm{CO}}_2$ 
% atoms, run:

tco2 = find(any(ATN.A(ico2,:), 1))
ATN.rxns(tco2)'
%% 
% i.e., three atom transitions in each of the reactions R3 and R4 involve atoms 
% in ${\textrm{CO}}_2$. To find atoms connected to ${\textrm{CO}}_2$ atoms via 
% these atom transitions, run:

cco2 = find(any(ATN.A(:, tco2) < 0,2));
ATN.mets(cco2)'
%% 
% i.e., ${\textrm{CO}}_2$ atoms are connected to atoms in the metabolites L-DOPA 
% (VMH ID: 34dhphe) and formate (VMH ID: for).
%% 
% *Step 2: Identify conserved moieties in DAS by graph theoretical analysis 
% of the atom transition network generated in Step 1.*

tic
[L,Lambda,moietyFormulas,moieties2mets,moieties2vectors,atoms2moieties] = ...
    identifyConservedMoieties(model, ATN);
t = toc;
fprintf('Computation time: %.1e s\n\n', t); % Print computation time
%% 
% This function outputs the moiety matrix (|L|), the moiety supergraph (|Lambda|), 
% the chemical formulas of moieties (|moietyFormulas|), and three vectors that 
% map between the various inputs and outputs. The 10×5 moiety matrix |L| has a 
% row for each metabolite and a column for each conserved moiety in DAS. Each 
% column is a moiety vector, with elements corresponding to the number of instances 
% of a conserved moiety in each metabolite. To find the number of instances of 
% moiety 2 in L-DOPA, run

iLDOPA = find(ismember(model.mets, '34dhphe[c]'))
full(L(iLDOPA, 2))
%% 
% i.e., L-DOPA contains one instance of moiety 2.
% 
% The 19×17 moiety supergraph (|Lambda|) contains the graphs of all seven conserved 
% moieties in DAS (Figure 7).
% 
% 
% 
% Figure 7: Graphs of the five conserved moieties in DAS. Each node represents 
% an instance of a conserved moiety in a particular metabolite. Each directed 
% edge represents conservation of a moiety between two metabolites. The chemical 
% formula of each moiety is given below its graph.
% 
% Each row of Lambda represents a single instance of a conserved moiety in a 
% particular metabolite. The vector moieties2vectors maps between the rows of 
% Lambda and the columns of L. To obtain the incidence matrix of a particular 
% moiety graph, e.g., λ2 in Figure 7, run

i2 = find(moieties2vectors == 2);
c2 = find(any(Lambda(i2, :)));
lambda2 = full(Lambda(i2, c2))
%% 
% The vector |moieties2mets| maps the rows of Lambda to metabolite indices in 
% the DAS reconstruction (|model|). To find metabolites containing instances of 
% moiety 2, run

m2 = moieties2mets(i2);
mets2 = model.mets(m2)'
%% 
% The chemical formula of moiety 2 is given by,

moietyFormulas{2}
%% 
% Finally, the vector atoms2moieties maps each atom in the atom transition network 
% for DAS to a particular instance of a conserved moiety. To find atoms in L-DOPA 
% that belong to moiety 2, run

find(ismember(atoms2moieties, i2) & ismember(ATN.mets, '34dhphe[c]'))'
%% 
% *Step 3: Classify moieties*

types = classifyMoieties(L, model.S)
%% 
% The internal moiety (λ3 in Figure 3) is conserved in both the open and closed 
% DAS network, whereas the transitive and integrative moieties are only conserved 
% in the closed network${\;}^6$.
%% Part 2: Effects of variable atom mappings between recurring metabolite pairs
% Here, we will again identify conserved moieties in DAS but with a slightly 
% different set of atom mappings (Figure 8). The different atom mappings gives 
% rise to a different atom transition network with a different set of conserved 
% moieties. In particular, it contains a single composite moiety, λ8 in Figure 
% 5, in place of the two moieties λ4 and λ5 in Figure 3. The composite moiety 
% is the result of variable atom mappings between the recurring metabolite pair 
% O2 and H2O in reactions R1 and R2.
% 
% 
% 
% Figure 8: (a) Oxygen atom transitions used in Part 1. Oxygen atom 1 in O2 
% maps to the oxygen atom in H2O in both R1 and R2. These atom transitions contain 
% two separate moieties, with two disconnected moiety graphs (λ4 and λ5 in Figure 
% 7), and two linearly independent moiety vectors (L(:,4) and L(:,5)). (b) Oxygen 
% atom transitions used in Part 2. A different oxygen atom maps from O2 to H2O 
% in R1 than in R2. These atom transitions contain only one composite moiety. 
% (c) The composite moiety graph arising from the oxygen atom transitions in (b).
%% 
% *Step 1: Identify conserved moieties with the alternative set of atom mappings.*

% Create an alternative MDL RXN file
R2rxn = regexp(fileread([outputDir filesep 'atomMapped' filesep 'R2.rxn']), '\n', 'split')';
R2rxn{2} = 'alternativeR2';
R2rxn{135}(62:63) = '18';
R2rxn{151}(62:63) = '19';
fid2 = fopen([outputDir filesep 'atomMapped' filesep 'alternativeR2.rxn'], 'w');
fprintf(fid2, '%s\n', R2rxn{:});
fclose(fid2);

% Create an alternative DAS model
alternativeModel = model;
alternativeModel.rxns{2} = 'alternativeR2';

% Identify conserved moieties
ATN = buildAtomTransitionNetwork(alternativeModel, atomMappedDir);
[L,Lambda,moietyFormulas,moieties2mets,moieties2vectors,atoms2moieties] = ...
    identifyConservedMoieties(alternativeModel, ATN);
%% 
% *Step 2: Decompose the composite moiety vector*
% 
% First, extract the internal stoichiometric matrix for DAS, by running:

rbool = ismember(alternativeModel.rxns, ATN.rxns);
mbool = any(alternativeModel.S(:,rbool), 2);
N = alternativeModel.S(mbool, rbool);
%% 
% To decompose the moiety matrix computed in Step 1, run:

try
    changeCobraSolver('gurobi6', 'milp');
end
D = decomposeMoietyVectors(L, N);
%% 
% Note that you can use any Mixed Integer Linear Programme (MILP) solver that 
% is supported by the COBRA toolbox. The decomposed moiety matrix D is identical 
% to the original moiety matrix computed in Part 1. Moiety vectors D(:,4) and 
% D(:,5) are the linearly independent components of the composite moiety vector 
% L(:,4) above.

full(D(:,[4 5])')
%% 
% One disadvantage of decomposing moiety vectors is that it is difficult to 
% keep track of which atoms belong to the decomposed moieties. We can, however, 
% estimate the chemical formulas of the decomposed moieties using the elemental 
% matrix for DAS. The elemental matrix is a numerical representation of the chemical 
% formulas of metabolites in DAS.

[E,elements] = constructElementalMatrix(alternativeModel.metFormulas,...
    alternativeModel.metCharges);
decomposedMoietyFormulas = estimateMoietyFormulas(D, E, elements);
decomposedMoietyFormulas([4 5])'
%% 
% i.e., each decomposed moiety contains an oxygen atom.
%% References
%% 
% # Haraldsdóttir, H.S., Thiele, I., Fleming, R.M. Comparative evaluation of 
% open source software for mapping between metabolite identifiers in metabolic 
% network reconstructions: application to Recon 2. _J. Cheminform_ 6(1), 2 (2014).
% # Elizabeth Brunk, et al. Recon 3D: A Three-Dimensional View of Human Metabolism 
% and Disease. Submited
% # Rahman, S.A., et al. Reaction Decoder Tool (RDT): extracting features from 
% chemical reactions.  _Bioinformatics_ 32(13), 2065–2066 (2016).
% # Preciat et al. Comparative evaluation of atom mapping algorithms for balanced 
% metabolic reactions: application to Recon 3D. _J Cheminform_, 9: 39 (2017).
% # Hulda S. Haraldsdóttir and Ronan M. T. Fleming. Identification of conserved 
% moieties in metabolic networks by graph theoretical analysis of atom transition 
% networks. _PLOS Comput. Biol_, 12(11) (2016).
% # Iman Famili and B. Ø. Palsson. The convex basis of the left null space of 
% the stoichiometric matrix leads to the definition of metabolically meaningful 
% pools. ‎_Biophys. J_, 85(1):16–26 (2003).
##### SOURCE END #####
-->
</div></body></html>